import "babel-polyfill";
import {map, reduce, range} from "../../../algorithm/iterator";
import Integer from "../../number/integer";
import {gcd as gcdInteger, invmod} from "../../number-theory/gcd";
import {generate} from "../../number-theory/prime";
import {sqrt} from "../../number-theory/nth-root";
import PolynomialZ, {CannotDivideError} from "../polynomial-Z";
import PolynomialZn from "../polynomial-Zn";
import {gcd as gcdPolynomial} from "../gcd";
import {factorUnique as factorZpUnique} from "./factor-Zp";
import LiftingTree from "./lifting-tree";

var squareFree = function*(polynomial) {
    var instance = PolynomialZ.instance;
    var f = polynomial.primitive();
    var h = PolynomialZ.zero;
    for (let degree = 0; ; ++degree) {
        let f1 = PolynomialZ::gcdPolynomial(f, f.deriv());
        let h1 = instance.divide(f, f1);
        let g = instance.divide(h, h1);
        if (degree > 0 && g.degree > 0) {
            yield [g, degree];
        }
        f = f1;
        h = h1;
        if (h.degree === 0) {
            break;
        }
    }
};

var selectPrime = function(polynomial) {
    for (let prime of generate()) {
        if (Integer.mod(polynomial.leadingCoeff, prime).sign) {
            let PolynomialType = PolynomialZn(prime);
            let p = new PolynomialType(polynomial);
            if (PolynomialType::gcdPolynomial(p, p.deriv()).degree === 0) {
                return prime;
            }
        }
    }
};

var norm2 = function() {
    return Integer.add(sqrt(this::map(Integer.square)::reduce(Integer.add)), Integer.one);
};

var landauMignotte = function() {
    return Integer.shl(Integer.multiply(this::norm2(), Integer.abs(this.leadingCoeff)), this.degree - 1);
};

var combination = function*(polynomial, fac, mod) {
    var check_and_divide = function(p, norm, degree) {
        if (Integer.compare(p::map(Integer.abs)::reduce(Integer.add), Integer.shl(norm, degree)) > 0) {
            return;
        }
        try {
            let [q, r] = polynomial.divmod(polynomial, p);
            if (r.degree === -Infinity) {
                return [p, q];
            }
        } catch (e) {
            if (!(e instanceof CannotDivideError)) {
                throw e;
            }
        }
    };
    var pre_check = function(norm, degree) {
        var coeff = reduce_half(Integer.mod(Integer.multiply(index::map(function(n) {
            var p = fac[n];
            return p.get(p.degree - 1);
        })::reduce(Integer.add), polynomial.leadingCoeff), mod));
        coeff = Integer.exactDivide(coeff, gcdInteger(coeff, polynomial.leadingCoeff));
        return Integer.compare(Integer.abs(coeff), Integer.pow(norm, degree)) < 0;
    };
    var product = function() {
        var p = new PolynomialZ(polynomial.leadingCoeff);
        for (let x of index) {
            p = new PolynomialZ(...p.multiply(p, fac[x])::map(n => Integer.mod(n, mod)));
        }
        return new PolynomialZ(...p::map(reduce_half)).primitive();
    };
    var combination_found = function() {
        norm = Math.min(norm, polynomial::norm2());
        var length = index.length;
        var fac2 = [];
        for (let i = 0; i < fac.length; ++i) {
            if (index.indexOf(i) === -1) {
                fac2.push(fac[i]);
            }
        }
        fac = fac2;
        if (index[0] + length - 1 >= fac.length) {
            index = [...range(length + 1)];
        } else {
            for (let i = 1; i < length; ++i) {
                index[i] = index[0] + i;
            }
        }
    };
    var next_combination = function(n) {
        var length = index.length;
        var m = length - 1;
        while (index[m] === fac.length - length + m) {
            --m;
        }
        if (m === -1) {
            index = [...range(length + 1)];
        } else {
            ++index[m];
            for (let i = m + 1; i < length; ++i) {
                index[i] = index[m] + i - m;
            }
        }
    };

    var half_mod = Integer.shr(mod, 1);
    var reduce_half = n => Integer.compare(n, half_mod) <= 0 ? n : Integer.subtract(n, mod);
    var norm = polynomial::norm2();
    var index = [0];
    while (index.length <= fac.length >>> 1) {
        let degree = index::map(x => fac[x].degree)::reduce((x, y) => x + y);
        let q = pre_check(norm, degree) && check_and_divide(product(), norm, degree);
        if (q) {
            yield q[0];
            polynomial = q[1];
            combination_found();
        } else {
            next_combination();
        }
    }
    yield polynomial;
};

export var factorUnique = function*(polynomial) {
    polynomial = PolynomialZ.instance.from(polynomial);
    var prime = selectPrime(polynomial);
    var PolynomialType = PolynomialZn(prime);
    var leading_coeff = polynomial.leadingCoeff;
    var p_poly = new PolynomialType(PolynomialZ.instance.multiply(polynomial, invmod(leading_coeff, prime)));
    var fac = [...PolynomialType::factorZpUnique(p_poly)];
    if (fac.length === 1) {
        yield polynomial;
        return;
    }
    var bound = polynomial::landauMignotte();
    var tree = new LiftingTree(polynomial, fac);
    var mod = prime;
    while (Integer.compare(mod, bound) < 0) {
        tree.lift();
        mod = Integer.square(mod);
    }
    fac = tree.liftResult::map(p => PolynomialZ.instance.from(p));
    mod = Integer.pow(prime, Math.floor((bound.bitLength - 1) / Integer.subtract(prime, Integer.one).bitLength));
    while (Integer.compare(mod, bound) < 0) {
        mod = Integer.multiply(mod, prime);
    }
    fac = [...tree.liftResult()::map(p => new PolynomialZ(...p::map(n => Integer.mod(n, mod))))];
    yield* combination(polynomial, fac, mod);
};

export default function* factor(polynomial) {
    polynomial = PolynomialZ.instance.from(polynomial);
    if (polynomial.degree === -Infinity) {
        yield [PolynomialZ.zero, 1];
        return;
    } else {
        let g = gcdInteger(...polynomial);
        if (polynomial.leadingCoeff.sign < 0) {
            g = Integer.negate(g);
        }
        if (!Integer.eq(g, Integer.one)) {
            yield [new PolynomialZ(g), 1];
        }
        for (let [p, degree] of squareFree(polynomial)) {
            for (let q of factorUnique(p)) {
                yield [q, degree];
            }
        }
    }
};
